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Modern society depends on increasingly interdepen- 
dent systems that are prone to widespread failure. Trans- 
portation, communication, power grids and other infras- 
tructures support one another and the world's inter- 
connected economies. Barrages of incidents large and 
small — downed power lines, grounded aircrafts, natural 
disasters and the like — cause avalanches of repercussions 
that cascade within and among these systems [T]. Al- 
though interdependence confers benefits, its effect on the 
risks of individual systems and on the collection of them 
remains poorly understood. 

Here we analyze how the interconnectivity (interde- 
pendence) between networks affects the sizes of their 
cascades of load shedding. For networks derived from 
interdependent power grids, we show that interdepen- 
dence can have a stable equilibrium. An isolated net- 
work suppresses its large cascades by connecting to other 
networks, but too many interconnections exacerbate its 
largest cascades — and those of the whole system. We de- 
velop techniques to estimate this optimal amount of in- 
terconnectivity, and we examine how differences among 
networks' capacity and load affect this equilibrium. Our 
framework advances the current mathematical tools for 
analyzing dynamics on interdependent (or modular) net- 
works, and it improves our understanding of systemic risk 
in coupled networks. 

In the basic process we consider, a system contains 
many elements that shed load to neighboring elements 
whenever they reach their capacity. This is captured 
by the classic sandpile model of Bak-Tang-Wiesenfeld, 
a paradigm for the power law statistics of cascades in 
many disciplines, from neuronal avalanches to financial 
instabilities to electrical blackouts [21]. In a basic for- 
mulation on a graph of nodes and edges, each node has 
a capacity for holding grains of sand (interpreted here 
as load or stress) . Grains of sand are dropped randomly 
on nodes, and whenever a node has more sand than its 
capacity, it topples and sheds all its sand to its neigh- 
bors, which may in turn have too much sand and topple, 
and so on. Thus dropping a grain of sand can cause an 
avalanche (cascade) of topplings. These avalanches, like 
blackouts in power grids |6J, occur in sizes characterized 



by a power law: they are often small but occasionally 
enormous. 

The Bak-Tang-Wiesenfeld model was originally formu- 
lated on a lattice. Given the relevance of networked sys- 
tems, the dynamics have recently been studied on iso- 
lated networks, but not yet on interdependent (or mod- 
ular) networks. Here we study it on two networks with 
sparse connections between them. Each network models 
an infrastructure (or a module of one), and the inter- 
connections between them model their interdependence. 
We explicitly study networks extracted from two inter- 
dependent power grids in the southeastern USA and an 
idealization of them that is more amenable to mathemat- 
ical study. In this idealization, each node is connected to 
a node in the other network with probability p (Fig. 1, 
inset). 
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FIG. Al: The chance that a network a coupled to another 
network h suffers a cascade larger than half its network (gold 
curve) has a stable minimum at a critical amount of inter- 
connectivity p* . Networks seeking to mitigate their largest 
cascades would prefer to build or demolish interconnections 
to operate at this critical point p*. The blue (red) curve is 
the chance that a cascade that begins in a (b) topples at least 
1000 nodes in a. Increasing interconnectivity only exacerbates 
the cascades inflicted from 6 to a (red), but interestingly it 
initially suppresses the local cascades in a. (From simulations 
on coupled random 3-regular graphs; the inset depicts a small 
example with 30 nodes per network and p = 0.1.) 
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Our main result is that interdependence can have a 
stable equilibrium (Fig. All. Some interconnectivity is 
beneficial to an individual network, for the other net- 
work acts as a reservoir for extra load. The gold curve 
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of Fig. XT shows that the chance of a large cascade in a 
network can be reduced by 70% by increasing the inter- 
connectivity p from 0.0005 to 0.075. Too much interde- 
pendence, however, becomes detrimental for two reasons. 
First, new interconnections open pathways for the neigh- 
boring network to inflict additional load. Second, each 
interconnection augments the system's capacity, making 
available more load that fuels even larger cascades in each 
network. As a result, the chance of a large cascade in an 
individual network eventually increases with interconnec- 
tivity p, so p* is a stable minimum. 

This second factor above — that new interconnections 
increase the networks' capacity for load — has global con- 
sequences. With more load available, larger cascades in 
the system as a whole become possible. Therefore net- 
works that interconnect to one another to mitigate their 
own cascades (Fig. |Al[ ) may unwittingly cause larger 
global cascades in the whole system. This is a warn- 
ing for the interconnections under construction among, 
for example, different power grids to accommodate long- 
distance trade and renewable sources of energy [TT] . 

The results in Fig. |A1| show that networks suppress- 
ing their largest cascades would seek interconnectivity p* . 
However, as shown in the the main article, building inter- 
connections to operate at p* increases the occurrence of 
small cascades. Conversely, networks can suppress their 
smallest cascades the most by seeking isolation, p = 0. 
But suppressing their smallest cascades exacerbates their 
largest ones (left side of Fig. All, just as extinguishing 



small forest fires can incite large ones and engineering 
power grids to suppress small blackouts can increase the 
risk of large ones [5] . 

Finally we determine how asymmetry among networks 
affects the optimal level of interconnectivity that each 
prefers. For instance, two interconnected power grids 
may differ in capacity, load, redundancies, demand, sus- 
ceptibility to line outages, and ages of infrastructure. We 
capture these differences with a parameter that controls 
the rates at which cascades begin in either network. We 
show that in any asymmetric situation the equilibrium 
will be frustrated, with only one of the grids able to 
achieve its optimal level of interconnectivity. 

Determining how interdependence affects the function- 
ing of networks is critical to understanding the infras- 
tructure so vital to modern society. Whereas others have 
recently shown that interdependence can lead to alarm- 
ingly catastrophic cascades of failed connectivity |20j . 
here we show that interdependence also provides ben- 
efits, and these benefits can balance the detriments at 
stable equilibria. We expect that this work will stim- 
ulate calculations of critical points in interconnectivity 
among networks subjected to other dynamics. As crit- 
ical infrastructures such as power grids, transportation, 
communication, banks and markets become increasingly 
interdependent, resolving the risks of large cascades and 
the incentives that shape them becomes ever more im- 
portant. 
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Suppressing cascades of load in interdependent networks 

Understanding how interdependence among systems affects cascading behaviors is increasingly important 
across many fields of science and engineering. Inspired by cascades of load shedding in coupled electric 
grids and other infrastructure, we study the Bak-Tang-Wiesenfeld sandpile model on modular random 
graphs and on graphs based on actual, interdependent power grids. Starting from two isolated networks, 
adding some connectivity between them is beneficial, for it suppresses the largest cascades in each system. 
Too much interconnectivity, however, becomes detrimental for two reasons. First, interconnections open 
pathways for neighboring networks to inflict large cascades. Second, as in real infrastructure, new inter- 
connections increase capacity and total possible load, which fuels even larger cascades. Using a multitype 
branching process and simulations we show these effects and estimate the optimal level of interconnectivity 
that balances their tradeoffs. Such equilibria could allow, for example, power grid owners to minimize the 
largest cascades in their grid. We also show that asymmetric capacity among interdependent networks 
affects the optimal connectivity that each prefers and may lead to an arms race for greater capacity. 
Our multitype branching process framework provides building blocks for better prediction of cascading 
processes on modular random graphs and on multi-type networks in general. 

I 



Networks that constitute our critical infrastructure in- 
creasingly depend on one another, which enables cascades 
of load, stress and failures [IHB]- The water network, for 
instance, turns turbines and cools nuclear reactors in the 
electrical grid, which powers the transportation and com- 
munications networks that underpin increasingly inter- 
dependent global economies. Barrages of disturbances at 
different scales — volcanic eruptions [7] , satellite malfunc- 
tions RP, earthquakes, tsunamis, wars [8] — trigger cas- 
cades of load shedding in interdependent transportation, 
communication and financial systems. Interdependence 
can also increase within a particular infrastructure. The 
electrical grid of the United States, for example, consists 
of over 3,200 independent grids with distinct ownership — 
some private, others public — and unique patterns of con- 
nectivity, capacities and redundancies [9]. To accommo- 
date rising demand for electricity, long distance trade of 
energy |lU], and new types of power sources jllj, inter- 
connections among grids bear ever more load |12) . and 
many new high-capacity transmission lines are planned 
to interconnect grids in the United States and in Eu- 
rope [13 . Figure [l] shows the new interconnections 
planned to transport wind power |llj . Though neces- 
sary, these interconnections affect systemic risk in ways 
not well understood, such as in the power grid, where the 
modular structure affects its large cascades. For exam- 
ple, the August 14, 2003 blackout, the largest in North 
American history, spread from a grid in Ohio to one in 
Michigan, then to grids in Ontario and New York before 
overwhelming the northeast [ini [H] . 

Researchers have begun to model cascades of load and 
failure within individual power grids using probabilistic 
models [5], linearized electric power dynamics [HI |T7] 
and game theory fl8^. The first models of interdepen- 
dent grids use simplified topologies and global coupling 
to find that interconnections affect critical points of cas- 
cades [19] , which suggests that they may affect the power 
law distributions of blackout size [6l [T2]- Models with 
interconnections among distinct infrastructure have fo- 
cused on the spread of topological failures, in which nodes 




FIG. 1: The power grid of the continental United States, illus- 
trating the three main regions or "interconnects" — Western, 
Eastern and Texas — and new lines (in orange) proposed by 
American Electric Power to transport wind power. Source: 
NPR [TT]. 



are recursively removed |20H22j , and not on the dynami- 
cal processes occurring on these networks. These models 
find that interdependence causes alarmingly catastrophic 
cascades of failed connectivity [20^22] . Yet as we show 
here, interdependence also provides benefits, and these 
benefits can balance the dangers at stable critical points. 

Here we develop a simple, dynamical model of load 
shedding on sparsely interconnected networks. We study 
Bak-Tang-Wiesenfeld (BTW) sandpile dynamics [23 Hi] 
on networks derived from real, interdependent power 
grids and on sparsely coupled, random regular graphs 
that approximate the real topologies. Sandpile dynam- 
ics are paradigms for the cascades of load, self-organized 
criticality and power law distributions of event sizes that 
pervade disciplines, from neuronal avalanches [25 -27 to 
cascades among banks [53] to earthquakes [53] , landslides 
[30] , forest fires [ST] [S^ , solar fiares [33] |34] , and electri- 
cal blackouts [6] . Sandpile cascades have been extensively 
studied on isolated networks [7HT^ [55] . On interdepen- 
dent (or modular) networks, more basic dynamical pro- 
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cesses have been studied [i^HiS] . but sandpile dynamics 
have not. 

We use a multitype branching process approximation 
and simulations to derive at a heuristic level how interde- 
pendence affects cascades of load. Isolated networks can 
mitigate their largest cascades by building interconnec- 
tions to other networks, as those networks provide reser- 
voirs to absorb excess load. Build too many interconnec- 
tions, however, and the largest cascades in an individual 
network increase in frequency for two reasons: neigh- 
boring networks inflict load more easily, and each added 
interconnection augments the system's overall capacity 
and load. These stabilizing and destabilizing effects bal- 
ance at a critical amount of interconnectivity, which we 
analyze for synthetic networks that approximate interde- 
pendent power grids. As a result of the additional load in- 
troduced by interconnections, the collection of networks, 
viewed as one system, suffers larger global cascades — a 
warning for the increasing interdependence among elec- 
trical grids (Fig. [l]), financial sectors and other infras- 
tructure pi¥n] . Finally we study the effects of capacity 
and load imbalance. Networks with larger total capac- 
ity inflict larger avalanches on smaller capacity networks, 
which suggests an arms race for greater capacity. The 
techniques developed here advance the theoretical ma- 
chinery for dynamical processes on multi-type networks 
as well as our heuristic understanding of how interde- 
pendence and incentives affect large cascades of load in 
infrastructure. 



I. THEORETICAL FORMULATION 
A. Sandpile dynamics 

Introduced by Bak, Tang and Wiesenfeld in 1987 and 
1988 [231 m], the sandpile model is a well-studied, styl- 
ized model of cascades that exhibits self-organized criti- 
cality, power laws and universality classes and that has 
spawned numerous related models with applications in 
many disciplines (e.g., [501 1311 1321 1311 ES])- In a basic 
formulation on an arbitrary graph of nodes and edges, 
one drops grains of sand (or "load") uniformly at ran- 
dom on the nodes, each of which has an innate threshold 
(or capacity). Whenever the load on a node exceeds its 
threshold, that node topples, meaning that it sheds (or 
moves) its sand to its neighbors. These neighbors may in 
turn become unstable and topple, which causes some of 
their neighbors to topple, and so on. In this way, drop- 
ping a grain of sand on the network can cause a cascade 
of load throughout the system — often small but occa- 
sionally large. The cascade finishes once no node's load 
exceeds its capacity, whereupon another grain of sand is 
dropped, and the process repeats. Probability measures 
of the size, area, and duration of avalanches typically 
follow power laws asymptotically in the limit of many 
avalanches [47] . 

Many classic versions of the sandpile model [531 ESI HZ] 



connect the nodes in a finite, two-dimensional lattice and 
assign all nodes threshold four, so that a toppled node 
sheds one sand grain to each of its four neighbors. The 
lattice has open boundaries, so that sand shed off the 
boundary is lost, which prevents inundation of sand. A 
few variants of the model on lattices can be solved exactly 
if the shedding rules have abelian symmetry I47j . 

More recently, sandpile models have been studied on 
isolated networks, including Erdos-Renyi graphs [TUII35] . 
scale-free graphs [ZHH], and graphs generated by the 
Watts-Strogatz model on one-dimensional [11] and on 
two-dimensional [T^ lattices. A natural choice for the 
capacities of nodes — which we use here — is their degree, 
so that toppled nodes shed one grain to each neighbor 
[II HHj . Other choices include identical [7| , uniformly dis- 
tributed from zero to the degree k [5], and k^^^ for some 
< ?7 < 1 [H1I3, but all such variants must choose ways to 
randomly shed to a fraction of neighbors. Shedding one 
grain to each neighbor is simpler and exhibits the rich- 
est behavior [7] . A natural analog of open boundaries on 
finite lattices is to delete grains of sand independently 
with a small probability / as they are shed. We choose 
the dissipation rate of sand f so that the largest cascades 
topple almost the entire network. 

The mean-field solution of sandpile cascades is charac- 
terized by an avalanche size distribution that asymptot- 
ically obeys a power law with exponent -3/2 and is quite 
robust to network structure. (For example, on scale-free 
random graphs, sandpile cascades deviate from mean- 
field behavior only if the degree distribution has a suffi- 
ciently heavy tail, with power law exponent 2 < 7 < 3 
[7].) Nevertheless, sparse connections among interdepen- 
dent networks divert and direct sandpile cascades in in- 
teresting, relevant ways. 



B. Topologies of interacting networks 

Here we focus on interdependent power grids and ide- 
alized models of them. We obtained topological data on 
two interdependent power grids — which we label c and 
d — from the US Federal Energy Regulation Commission 
(FERC) [49]. (All data shown here is sanitized to omit 
sensitive information.) Owned by different utilities but 
connected to one another in the southeastern USA. [55] 
power grids c and d have similar size (439 and 504 buses) 
but rather different average internal degrees (2.40 and 
2.91, respectively). The grids are sparsely interconnected 
by just eight edges, making the average external node 
degrees 0.021 and 0.018, respectively. More information 
on c, d is in Table SI of the SI. As in other studies, we 
find that these grids have narrow degree distributions 
[31 [Sni [5T] and are sparsely interconnected to one an- 
other [12 . 

To construct idealized versions of the real grids, con- 
sider two networks labeled a and b. Due to the narrow 
degree distribution of the real grids, we let network a be 
a random Za-regular graph (where each node has degree 



5 



Za) and network 6 be a random 2:(,-regular graph. These 
two are then sparsely interconnected as defined below. To 
define this system of coupled networks more formally, we 
adopt the multitype network formalism of [iSJ [S3] . Each 
network a, h has its own degree distribution, Paikaa, kab) 
and Pb{kba,kbb), where, for example, Pa{kaa-,kab) is the 
fraction of nodes in network a with fc^a neighbors in a and 
kab in b. We generate realizations of multitype networks 
with these degree distributions using a simple generaliza- 
tion of the configuration model: all nodes repeatedly and 
synchronously draw degree vectors (fcoai^ob) from their 
degree distribution po (where o G {a, &}), until the totals 
of the internal degrees kaa, kbb are both even numbers and 
the totals of the external degrees kab, kba are equal. [69] 

We interconnect the random 2{,-regular graphs by 
Bernoulli- distributed coupling: each node receives one 
external "edge stub" with probability p and none with 
probability 1 — p. Hence the degree distributions are 

Pa{ZaA) = P,Pa{Za,Q) = 1 - P, and Pbil, Zb) = 

p,Pb{0,Zb) = I — p. We denote this class of interacting 
networks by the shorthand R{za)-B{p)-R{zb); we illus- 
trate a smaU example of i?(3)-B(0.1)-i?(4) in Fig. [2j 




FIG. 2: Random 3- and 4-regular graphs connected by 
Bernoulli-distributed coupling with interconnectivity param- 
eter p = 0.1 {R{3)-B{0.l)-R{4)). We also illustrate the shed- 
ding branch distribution qab(rba,bb), the chance that an ab- 
shedding causes rba,rbb many bo-, fefe-sheddings at the next 
time step. Note these random graphs are small and become 
tree-like when they are large 1000 nodes). 



C. Measures of avalanche size 

We are most interested in the avalanche size distribu- 
tions Sa{ta-,tb),Sb{ta-,tb), whcre, for example, Sa{ta,tb) is 
the chance that an avalanche begun in network a (indi- 
cated by the subscript on Sa) causes ta, tb many topplings 
in networks a, b, respectively. These distributions count 
the first toppling event, and they are defined asymptoti- 
cally in that Sa, Sb are frequencies of avalanche sizes after 
the networks have undergone many cascades. To study 
Sa and Sbi we simulate sandpile avalanches and approxi- 
mate them using a multitype branching process. 



D. Multitype branching process approximation 

In these next two sections, we present an overview 
of our mathematical formulation, with details left to the 
Materials and Methods. We develop a branching pro- 
cess approximation that elucidates how sandpile cascades 
spread in interconnected networks, advances theoretical 
tools for cascades on multitype networks, justifies using 
this model as an idealization of real infrastructure like 
power grids, and establishes an open and relevant math- 
ematical challenge. However, readers more interested in 
the applications of the model may wish to skip to the 
Results section. 

Sandpile cascades on networks can be approximated 
by a branching process provided that the network is lo- 
cally tree-like (i.e., has few short cycles), so that branches 
of a nascent cascade grow approximately independently. 
The interacting networks R{za)-B{p)-R{zb) are tree-like 
provided they are sparse and large enough (with at least 
several hundred nodes), since the edges are wired uni- 
formly at random. Power grids are approximately tree- 
like: the clustering coefficient of power grids c and d, 
for example, is C « 0.05, an order of magnitude larger 
than an Erdds-Renyi random graph with equally many 
nodes and edges, but still quite small. Although tree- 
based approximations of other dynamical processes work 
surprisingly well even on non-tree-like graphs, power grid 
topologies were found to be among the most difficult to 
predict with tree-based theories [5]. Here we find that 
analytic, tree-based approximations of sandpile dynam- 
ics agree well with simulations even on the real power 
grid topologies (Figs, [s) S7, S8). 

Cascades on interacting networks require a multitype 
branching process, in which a tree grows according to 
probability distributions of the number of events of var- 
ious types generated from seed events. We consider two 
basic event types, a-topplings and 6-topplings — i.e., top- 
pling events in networks a and in h. These simplify the 
underlying branching process of sheddings, or grains of 
sand shed from one network to another, of which there 
are four types: aa-, ab-, ba- and 56-sheddings. (Note that 
there is no distinction between topplings and sheddings 
on one, isolated network, because sand can only be shed 
from, say, a to a.) 

A key property of sandpile dynamics on networks, 
which enables the branching process calculations, is that 
in simulations the amount of sand on a node is asymp- 
totically uniformly distributed from zero to one less its 
degree (i.e., there is no typical amount of sand on a 
node) jlSl [SB] . Hence the chance that a grain arriving 
at a node with degree k topples it equals the chance that 
the node had k — 1 grains of sand, which is 1/k. So 
sandpile cascades are approximated by what we call 1/k- 
percolation: the cascade spreads from node u to node v 
with probability inversely proportional to the degree of 
V. This suggests a direct interpretation for infrastruc- 
ture: important nodes have k times more connectivity 
than unimportant (degree-1) nodes, so they are k times 
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less likely to fail (they are presumably reinforced by en- 
gineers). But when important nodes do fail, they cause 
k times more repercussions (shedded grains of sand) . We 
found some evidence for this in the power flowing through 
buses (nodes) in power grids: each additional degree cor- 
relates with an additional 124 MVA of apparent power 
flowing through it (i?^ = 0.30; see Fig. S6 of the SI). 

The details of the branching process analysis extend 
the standard techniques as presented in the Materials and 
Methods. We give here only the crux of the derivation. 
Suppose a grain of sand is shed from network o € {a, b} 
to network d S {a,b} ('o' for "origin network", 'd' for 
"destination network"). What is the chance that this 
grain shed from o to d (an od-shedding) causes r^a and 
rj^h many grains to be shed from network d to a and from 
d to &, respectively, at the next time step? This probabil- 
ity distribution, denoted qod{rda,rdb), is the branch (or 
children) distribution of the branching process for shed- 
dings. Figure [2] illustrates qab as an example. Neglect- 
ing degree-degree correlations (the subject of so-called 
P{k,k') theory [5]), a grain shed from network o to d 
arrives at an edge stub chosen uniformly at random, so 
it arrives at a node with degree Pd(»'rfa, ''df)) with proba- 
bility proportional to r^o, since that node has rdo many 
edges pointing to network o. Using this and the chance 
of toppling found above (1 /total degree), we approximate 
asymptotically that 



qodirda,rdb) 



rdoPd{rda,rdb) 



1 



(fc, 



do) 



Tda + Tdb 



(1) 



for Tda + Tdb > 0, where (kdo) is the expected number of 
edges from d to o, J2k^^,k^^ kdoPdikda, kdb)- To normaUze 
Qod, set 



9od(0,0) := 1 



E 

rda+rdb>0 



qod{rda,rdb), 



(2) 



which is the probability that the destination node does 
not topple (i.e., that it has fewer grains than one less its 
total degree). 

Note that for an individual, isolated network the anal- 
ogous branch distribution q[k) simplifies considerably: in 
the equivalent of Eq. ([ij there is a cancelation of k in the 
numerator of with 1/k on the right [7-9]. Thus the ex- 
pected number of children events {q{k)) — J2k ^^717^ i — 
1. Each seed event gives rise to one child on average, 
which then gives rise to one child on average, etc., which 
is called a "critical" branching process. (If less than one 
child on average, the branching process dies out; if more 
than one it may continue indefinitely.) The branching 
process approximations of sandpile cascades on the in- 
teracting networks studied here — coupled random reg- 
ular graphs a, b and power grids c, d — are also critical, 
because the principle eigenvalue of the matrix of first 
moments of the branch distributions is one. 

The branching process of sheddings is high dimen- 
sional, with four types aa, ab, ba, bb recording origin and 
destination networks. Transforming the shedding branch 



distributions qod to the toppling branch distributions 
Ua,Ub is easy; the key is that a node topples if and only 
if it sheds at least one grain of sand (for the details, see 
the Materials and Methods). This also halves the dimen- 
sions of the branching process of topplings, simplifying 
calculations. 



E. Self-consistency equations 

We analyze implicit equations for the avalanche size 
distributions Sa,Sb using generating functions [57^ De- 
note the generating functions associated to the toppling 
branch distributions u^, Ub and the avalanche size distri- 
butions Sq, Sb by capital letters lA and 5; for example. 



Ua{Ta,Tb) 



E 



Ua{ta,tb)Tl''Tl'' for Tda,Tdb G 



The theory of multitype branching processes [5 8) implies 

the self- consistency equations 

Sa ~ TaUa{Sa,Sb), Sb = TbUb{Sa, Sb) , (3) 

where each S is evaluated at (Ta,Tf,). In words, the left- 
hand equation in (|3| says that to obtain the distribution 
of the sizes of cascades begun in a, the cascade begins 
with an a-toppling (hence the Ta out front), which causes 
at the next time step a number of a- and 6-topplings 
distributed according to Ua, and these topplings in turn 
cause numbers of a- and fe-topplings distributed according 
to Sa and Sb. 

We wish to solve Eqs. (|3| for Sa and Sb, because their 
coefficients are the avalanche size distributions Sa,Sb of 
interest. In practice, however, these implicit equations 
are transcendental and difficult to invert. Instead, we 
solve Eqs. ([3]) with computer algebra systems using three 
methods — iteration, Cauchy integral formula, and multi- 
dimensional Lagrange inversion [13^ — to compute exactly 
hundreds of coefficients; for details, see the Materials and 
Methods. Figure [3] shows good agreement between sim- 
ulation of sandpile cascades on power grids c, d and the 
branching process approximation (obtained by iterating 
Eqs. ([3| seven times starting from Sa = Sb = I, with 
branch distributions calculated from the empirical degree 
distributions of c, d). For more details on the agreement, 
including how the degrees of nodes with external links ac- 
count for the characteristic "blips" in the avalanche size 
distributions of the power grids, see Figs. S7 and S8 of 
the SI. 

These numerical methods are computationally fea- 
sible for the probabilities of the smallest avalanches, 
but we are most interested in the probabilities of the 
largest avalanches — that is, in the asymptotic behav- 
iors of Sa(ta,tb), Sb(ta,tb) as ta,tb — >■ oo. Unfortunately 
the technique used for an isolated network — an expan- 
sion at a singularity of U — fails for sandpile cascades 
on Bernoulli-coupled random regular graphs and for the 
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FIG. 3: The multitype branching process (red curves) approx- 
imates simulations of sandpile cascades on the power grids 
c, d (blue curves) surprisingly well, given that power grids 
are among the most difhcult network topologies on which to 
predict other dynamics [5]. We plot the four marginalized 
avalanche size distributions in order to view one-dimensional 
curves (and here we label power grids c, d as a, b, respectively). 



power grids, because their generating functions have sin- 
gularities at infinity and none in the finite plane (see 
Materials and Methods). Generalizing these asymptotic 
techniques to "multitype cascades" with singularities at 
infinity poses an outstanding mathematical challenge. 
Nevertheless, three tactics — simulations, computer cal- 
culations of coefficients of Sa,Sb, and analytical calcula- 
tions of the first moments of the branch and avalanche 
size distributions — suffice to obtain interesting conclu- 
sions about the effect of interdependence on critical cas- 
cades of load, as discussed next. 



II. RESULTS 
A. Locally stabilizing effect of interconnections 

We first answer the question, would an isolated net- 
work suppress its largest cascades of load by connect- 
ing to another network? For coupled random regular 
graphs R{za)-B{p)-R{zh) , yes: increasing interconnectiv- 
ity p suppresses an individual network's largest cascades, 
but only up to a critical point p* (Fig. |4]). 

First we introduce notation. For a cascade that begins 
in network a, the random variables Taa,Tab are the sizes 
of the "local" and "inflicted" cascades: the number of 
topplings in a and in b, respectively. For example, a cas- 
cade that begins in a and that topples 10 a-nodes and 5 
6-nodes corresponds to Taa = 10, Tat = 5. We denote Tq 
to the be random variable for the size of a cascade in net- 
work a, without distinguishing where the cascade begins. 
(We define T(,a,Tbb,Tb analogously.) Dropping sand uni- 
formly at random on two networks of equal size means 
that avalanches begin with equal probability in either 

network, so Pr(Ta = ta) = J2tS^a{ta,tb) + Sb{ta,tb))/2. 

In Fig. |4] we plot the probability of observing a large 



avalanche in a (that topples at least half of all its nodes) 
as a function of interconnectivity p, as measured in nu- 
merical simulations on the R{3)-B{p)-R{3) topology. We 
distinguish between those avalanches that begin in a 
(blue "local cascades"), begin in b (red "inflicted cas- 
cades"), or in either network (gold). With increasing 
interconnectivity p, large inflicted cascades from a to 6 
(red curve) increase in frequency due in large part to 
the greater ease of avalanches traversing the intercon- 
nections between networks. More interesting is that in- 
creasing interconnectivity suppresses large local cascades 
(blue curve) for small p, but amplifies them for large 
p. The 80% drop in Pr(Taa > 1000) and 70% drop in 
PiiTa > 1000) from p ^ 0.001 to p* w 0.075 ± 0.01 
are the locally stabilizing effects of coupling networks. 
The left inset to Fig. [4] is the rank-size plot showing 
the sizes of the largest avalanches and their decrease 
with initially increasing p, and the same holds for sim- 
ulations on the power grids c and d (right inset). [70] 
The curve Pr(rci > C) and the location of its criti- 
cal point p* in Fig. [2] is robust to changing the cutoff 
C G [400,1500]. Thus a network such as a seeking to 
minimize its largest cascades would seek interconnectiv- 
ity that minimizes Pr(Ta > C), which we estimate to be 
p* w 0.075 ±0.01 for R{3)-B{p)-R{3) with 2 x 10^ nodes 
per network. 
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FIG. 4: Interconnectivity is locally stabilizing, but only up 
to a critical point. The main plot, the results of simula- 
tions on R(3)-B{p)-R{3) (2 x 10^ grains, / = .01, 2 x 10^ 
nodes/network), shows that large local cascades decrease and 
increase with p, while large inflicted cascades only become 
more likely. Their average (gold curve) has a stable mini- 
mum at p* ~ 0.075 ± 0.01. This curve and its critical point 
is stable to cutoffs c different from 1000 (400 < x < 1500). 
Inset figures: rank-size plots on log-log scales of the largest 
cascades in network a (left) for p = 10""^, 10~^, 10~^ and in 
power grid d connected to c by 0, 8 or 16 edges. Both cases 
are in the regime of sparse interconnectivity (left-hand side 
of main plot), as greater interconnectivity suppresses large 
cascades in an individual network. 



This central result appears to be generic: changing the 
system size, internal degrees, or type of degree distribu- 
tion (so long as it remains narrow) may slightly change 
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p* but not the qualitative shape of Fig. |4] (see SI Figs. 
SI, S2). Furthermore, this effect of optimal connectiv- 
ity is unique to interconnected networks: adding edges 
to a single, isolated network does not reduce the chance 
of a large cascade (and hence produce a minimum like 
in Fig. k|). (This is expected since 
narrow degree distributions [^.) Whereas adding links 
within a network can only increase its vulnerability to 
large cascades, adding links to another network can re- 
duce it. Note that we cannot derive analytical results 
for Fig. |4] because the standard techniques for single net- 
works fail for multitype generating functions with singu- 
larities at infinity, and inverting Eq. ([3| numerically is 
practical only for exactly computing the probabilities of 
small cascades (Tj, < 50) and not large ones (T^ > 10^) 
(see the Materials and Methods). These pose open and 
relevant mathematical challenges for future work. 

Intuitively, adding connections between networks di- 
verts load, and that diverted load tends to be absorbed 
by the neighboring network rather than amplified and re- 
turned, as most cascades in isolated networks are small. 
One way to see the diversion of load is in the first mo- 
ments of the toppling branch distributions UajUi,- For 
R{za)-B{p)-R{zb), the average numbers of topplings at 
the next time step in the same and neighboring networks, 
respectively, decrease and increase with the interconnec- 

tivity p as {Ua)a = 1 - P/(l + ^a), {Ua)b = P/(l + ^b), 

where, for example, {ua)a = J2t^,u ^aUaita^h)- 

However, introducing too many interconnections is 
detrimental as shown in Fig. [4j Interconnections let 
diverted load more easily return and with catastrophic 
effect. In addition, each interconnection augments the 
networks' capacity and hence average load, so that large 
avalanches increase in frequency in individual networks 
and in the collection of networks, as discussed next. 



B. Globally destabilizing effect of interconnections 

Adding interconnections amplifies large global cas- 
cades. That is, the largest avalanches in the collection 
of networks — viewed as one system with just one type 
of node — increase in size with interconnectivity. Here we 
are interested in the total avalanche size distribution s(t), 
the chance of t topplings overall in a cascade. Figure [5] 
shows the extension of the right-hand tail of s{t) in sim- 
ulations on R(S)-B{p)-R[3) with increasing interconnec- 
tivity p. The rank-size plot inset shows more explicitly 
that the largest avalanches increase with p. Similar re- 
sults occur in simulations on power grids c, d (see page 2 
of the SI). 

What amplifies the global cascades most significantly 
is the increase in total capacity (and hence average load 
available for cascades) and not the increased interdepen- 
dence between the networks. (Recall that capacities of 
nodes are their degrees, so introducing new edges be- 
tween networks increments those nodes' degrees and ca- 
pacities.) To see this effect on coupled random regular 




FIG. 5: Increasing the interconnectivity p between two ran- 
dom 3-regular graphs extends the tail of the total avalanche 
size distribution s{t), which does not distinguish whether 
nodes are in network a or 6 but considers them as one net- 
work. The inset shows a rank-size plot on log-log scales of 
the number of topplings t in the largest 10^ avalanches (with 
2 X 10^ grains of sand dropped) , confirming that adding more 
edges between random 3-regular graphs enlarges the largest 
global cascades by an amount on the order of the additional 
number of interconnections. 



graphs, we perform the following rewiring experiment. 
Beginning with two isolated random regular graphs, each 
node changes one of its internal edge stubs to be exter- 
nal with probability p. The degree distributions become, 
for example, Paiza - 1,1) = p,Pa{za,0) = I - p, which 
we call "Correlated-Bernoulli coupling" because the in- 
ternal and external degrees are not independent. Figure 
S4 of the SI shows that the largest global avalanches are 
not significantly enlarged with increasing "rewired inter- 
connectivity" p for random 3-regular graphs with such 
coupling. Furthermore, the enlargement of the largest 
cascades observed in the rank-size plot in the inset of 
Fig. [5] is on the order of the extra average load result- 
ing from the additional interconnectivity (and the same 
holds for simulations on the power grids). 

The amplification of global avalanches, though rel- 
atively small, is relevant for infrastructure: addi- 
tional capacity and demand often accompany — and even 
motivate — the construction of additional interconnec- 
tions |lllll3j . Furthermore, in reality it is more common 
to augment old infrastructure with new interconnections 
as in Fig. [l] (Bernoulli coupling) rather than to delete an 
existing internal connection and rewire it to span across 
networks (Correlated-Bernoulli coupling). Thus, build- 
ing new interconnections augments the entire system's 
capacity, and hence average load, and hence largest cas- 
cades. 
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C. Interconnect ivity that mitigates cascades of 
different sizes 

Figure |4] shows that networks seeking to mitigate 
their large avalanches seek optimal interconnectivity 
p. Networks mitigating their small or intermediate 
cascades would seek different optimal interconnectivity 
p, as shown in Fig. [6] (the results of simulations on 
R{3)-B{p)-R{3)). Figure [6j^ shows that the probabil- 
ity of a small cascade in network a {1 < Ta < 51) 
increases monotonically with p, so that networks miti- 
gating the smallest cascades seek isolation, p — 0. By 
contrast, the chance of a cascade of intermediate size, 
100 < Ta < 150 (Fig. [6|3), has a local maximum at 
p* « 0.05, so networks mitigating intermediate cascades 
would demolish all interconnections {p = 0) or build as 
many as possible {p — 1). By contrast, the largest cas- 
cades 400 <Ta < 1500 (Figs, [ep, [6p) occur with min- 
imal probability at p* « 0.075 ± 0.01. For more plots 
showing the change in concavity and the robustness of 
the stable critical point p* for large cascades, see Fig. 
S3. 



Pr(l<r„<51) Pr(100 < r,, < 150) 




FIG. 6: Networks mitigating the smallest cascades 1 < Ta < 
51 (A) seek isolation p = 0, while networks suppressing inter- 
mediate cascades 100 < Ta < 150 (B) seek isolation p = or 
strong coupling p = 1, depending on the initial interconnec- 
tivity p in relation to the unstable critical point p* ~ 0.05. 
But networks like power grids mitigating large cascades (C, 
D) would seek interconnectivity at the stable equilibrium 
p* ~ 0.075 ± 0.01. The bottom figures and the location of 
p* are robust to changes in the window i < Ta < £ -|- 50 for 
all 400 < ^ < 1500. 



Other models of cascades in power grids conclude that 
upgrading and repairing the system to mitigate the small- 
est blackouts may increase the risk of the largest black- 
outs 1^. Similarly, extinguishing small forest fires, a com- 
mon policy in the 20th century, increases forest vegeta- 
tion and thus the risk of large forest fires — a phenomenon 
known as the "Yellowstone effect" [32]. The results here 
augment these conclusions to include interconnectivity. 
Networks building interconnections in order to suppress 



their large cascades cause larger global cascades (by an 
amount on the order of the additional capacity). Net- 
works suppressing their small or intermediate cascades 
may seek isolation (p = 0), which amplifies their large 
cascades, or (for intermediate cascades) strong coupling 
{p — 1), which amplifies both large local and global cas- 
cades. 



D. Capacity disparity 

Two networks that are interdependent are rarely iden- 
tical, as are the R{3)-B{p)-R{3) topologies studied thus 
far, so we determine the effect of capacity disparity on 
cascades. As the capacities of nodes are their degrees, we 
study R{za)-B{p)-R{zb) with different internal degrees, 
Za 7^ We find that interdependence is more catas- 
trophic for smaller capacity networks, in that they suffer 
comparatively larger inflicted cascades. They still prefer 
some interconnectivity, but less than the higher capacity 
network. 

Using the theoretical branching process approxima- 
tion, we compute how much larger inflicted cascades 
are from high- to low-capacity networks. Differentiating 
Eqs. ([3| with respect to Ta,Ti, and setting = = I 
yields four linear equations for the first moments of 
the avalanche size distributions Sa, Sb in terms of the 
first moments of the branch distributions Ua,Ub. For 
R{za)-B{p)-R{zb), the four first moments of Sa,S{, are 
all infinite, as expected, because in isolation these net- 
works' avalanche size distributions are power laws with 
exponent —3/2 [THS]. Nevertheless, we can compare the 
rates at which the average infiicted cascade sizes diverge 
by computing their ratio 

{Sa)b ^ l + Zg 

{Sb)a l+Zb' ^' 

where, e.g., {sa)a = Y^t^^t^taSaita.tb) (see SI for the 
derivation). Thus Zb > Za implies that the inflicted cas- 
cades from 6 to a are larger on average than those from a 
to b. Fig. S9 of the SI shows qualitative agreement with 
simulations. 

As a result of Eq. Q , the low-capacity network prefers 
less interconnectivity than the high-capacity network. In 
simulations of i?(3)-i3(p)-i?(4), for instance, low capacity 
network a prefers p* ~ 0.05, whereas high-capacity b 
prefers pi w 0.3. For systems like power grids seeking to 
mitigate their cascades of load, these results suggest an 
arms race for greater capacity to fortify against cascades 
inflicted by neighboring networks. 



E. Incentives and equilibria in power grids 

Since different networks have unique susceptibilities to 
cascades (due to capacity disparity, for example), equi- 
libria among real networks are more nuanced than on 
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identical random graphs. Next we explore how the level 
of interconnectivity and load disparity affect sandpile 
cascades on the power grids c and d. (Although sand- 
pile dynamics do not obey Ohm's and Kirchoff's laws 
nor the flow of load from sources to sinks, as in phys- 
ical power flow models (e.g., [El E]), they do closely 
resemble some engineers' models of blackouts, and black- 
out data show evidence of criticality and power laws {&,.) 
To interpret results, we suppose that the owners of the 
power grids c, d are rational, in that they wish to miti- 
gate their largest cascades but care little about cascades 
overwhelming neighboring grids. 

To capture different amounts of demand, numbers of 
redundancies, ages of infrastructure, susceptibility to sag- 
ging power lines [111 [T7] , and other factors that affect the 
rate at which cascades of load shedding and failures begin 
in each network, we introduce a load disparity parameter 
r as follows. Each node in c is r times more likely than 
a node in d to receive a new grain of sand. Increasing r 
intensifies the load on grid c, the rate at which cascades 
begin there, and the sizes of the largest inflicted cascades 
from c to d. The larger r is, the more volatile power grid 
c becomes. 

Given the spatial structure of the power grid networks, 
there is no principled way to add arbitrarily many inter- 
connections between them. However, three different lev- 
els of interconnectivity are natural: (1) delete the eight 
interconnections so that c and d are isolated, (2) leave 
the eight original interconnections, and (3) add eight ad- 
ditional interconnections in a way that mirrors the em- 
pirical degree distribution. 




Rank 



FIG. 7: The critical load disparity at which inflicted cascades 
in d become equally large as its local cascades is r* ~ 15 
(for 16 interconnections ("bridges")). (For 8 interconnections, 
r* > 20; see Fig. SI SI. 4 1. Here we show a rank-size plot in 
log-log scales of the largest 10'' avalanches in power grid d, 
distinguishing whether they begin in c ("inflicted cascades") 
or in d ("local cascades"), for 8 and 16 interconnections (solid, 
dashed curves), for r = 1 (main), r — 15 (inset), in simula- 
tions with dissipation of sand / = 0.05, 10® grains dropped 
(after 10^ transients). 



Figure [7] shows that there are two ways to amplify the 



largest inflicted cascades in d. The first is to increase the 
number of interconnections (compare red to green-dashed 
curve). The second is to increase r (compare distance 
from green to gold curve in the main and in the inset). 
At a critical r* , the largest inflicted cascades in d that 
begin in c (red and green curves) are equally large as the 
largest local cascades in d that begin in d (blue and gold 
curves). For 16 interconnections, we estimate r* w 15 
(inset of Fig.j?]), and the inflicted cascades are larger and 
smaller, respectively, for r = 10, 20 (Fig. S5) — indicating 
that inflicted cascades begin to dominate local cascades 
at r* « 15. The actual load disparity between power 
grids c and d is r « 0.7, which we estimate by computing 
the average power incident per node in simulation output 
from FERC [49 . (There are, however, interdependent 
power grids in the southeastern USA with r > 15.) Since 
r — 0.7, the load is greater on grid d, so grid d prefers 
more interconnections and c prefers fewer than if r were 
1. Consequently, any equilibrium between the two grids 
is frustrated (or semi-stable): only one grid can achieve 
its desired interconnectivity. 



III. DISCUSSION 

We have presented a comprehensive study of sandpile 
cascades on interacting networks to obtain a deeper un- 
derstanding of cascades of load on interdependent sys- 
tems, showing both the benefits and dangers of inter- 
dependence. We combine a mathematical framework 
for multitype networks |451 153) with models of sand- 
piles on isolated networks [7HT^ to derive a multi- 
type branching process approximation of cascades of load 
between simple interacting networks and between real 
power grids. We show that some interdependence is ben- 
eficial to a network, for it mitigates its largest avalanches 
by diverting load to neighboring networks. But too much 
interconnectivity opens new avenues for infiicted load and 
adds capacity that fuels even larger cascades. The bene- 
fits and detriments in mitigating large avalanches balance 
at a critical amount of interconnectivity, which is a sta- 
ble equilibrium for coupled networks with similar capaci- 
ties. For coupled networks with asymmetric capacity, the 
equilibria will be frustrated, in that the networks prefer 
different optimal levels of interconnectivity. 

We also show that tuning interconnectivity to sup- 
press cascades of a certain range of sizes amplifies 
the occurrence of cascades in other ranges. Thus a 
network mitigating its small avalanches amplifies its 
large ones (Fig. |4]), and networks suppressing their own 
large avalanches amplify large ones in the whole system 
(Fig. [5]). Similarly it has been found that mitigating 
small electrical blackouts and small forest fires appears 
to increase the risk of large ones [51 [32]. Furthermore, 
the amplification of global cascades due to the increase 
in capacity (Fig.jsJ) is a warning for new interconnections 
among power grids (Fig. [T]) [TI] [T3] , financial sectors and 
other infrastructure. 
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These findings suggest economic and game-theoretic 
imphcations for infrastructure. (Note that here we con- 
sider the strategic interactions of entire networks rather 
than of individual nodes in one network, which is more 
standard, e.g., ,6QJ.) For example, a power grid owner 
wishing to mitigate the largest cascades in her grid would 
desire some interconnections to other grids, but not too 
many. However, what benefits individual grids can harm 
society: grids building interconnections amplify global 
cascades in the entire system. More detailed models that 
combine results like these with economic and physical 
considerations of electrical grids and with costs of build- 
ing connections may provide more realistic estimates of 
optimal levels of interconnectivity. Our framework — 
which models a dynamical process on stable, underlying 
network topologies — could also be combined with models 
of topological failure in interdependent networks [20H22] . 
Those studies conclude that systemic risk of connectiv- 
ity failure increases monotonically with interdependence 
("dependency links"). Whether suppressing cascades of 
load or of connectivity failures is more important might 
suggest whether to interconnect some {p* > 0) [Fig. [4] or 
none {p = 0) 20-22], respectively. Our results are con- 
sistent with recent work showing that an intermediate 
amount of connectivity minimizes risk of systemic de- 
fault in credit networks [HI], in contrast to Refs. 
and the more traditional view that risk increases mono- 
tonically with connectivity in credit networks (e.g., |62)). 

This work also advances our mathematical understand- 
ing of dynamical processes on multitype networks. Since 
networks with one type of node and edge are impover- 
ished views of reality, researchers have begun to study 
dynamical processes on multitype networks, such as on 
modular graphs [i2U45j . Here we derive a branching 
process approximation of sandpile cascades on multitype 
networks starting from the degree distributions, and we 
discuss open problems in solving for the asymptotic be- 
havior of the generating functions' coefficients, which 
elude current techniques for isolated networks. We ex- 
pect that the computational techniques used here to solve 
multidimensional generating function equations, such as 
multidimensional Lagrange inversion |13j . will find other 
uses in percolation and cascades in multitype networks. 
Finally, in the Appendix we derive the effective degree 
distributions in multitype networks, which expands the 
admissible degree distributions that others have consid- 
ered. The machinery we develop considers just two in- 
teracting networks, a and 6, or equivalently one network 
with two types of nodes. However, this extends to finitely 
many types, which may be useful for distinguishing types 
of nodes — such as buses, transformers and generators 
in electrical grids — or for capturing geographic informa- 
tion in a low-dimensional way without storing explicit 
locations — such as buses in the interiors of power grids 
and along boundaries between them. 

Here we have focused on mitigating large avalanches in 
modular networks, but other applications may prefer to 
amplify large cascades, such as the adoption of products 



in modular social networks |63j or propagating response- 
driven surveys across bottlenecks between social groups 
|64) . Cascades in social networks like these may require 
networks with triangles or other subgraphs added [531 1551 
I66j : inverting the resulting multidimensional generating 
function equations for dynamics on these networks would 
require similar multitype techniques. 

We expect that this work will stimulate calculations of 
critical points in interconnectivity among networks sub- 
jected to other dynamics, such as linearized power flow 
equations in electrical grids |16l 117] and other domain- 
specific models. As critical infrastructures such as power 
grids, transportation, communication, banks and mar- 
kets become increasingly interdependent, resolving the 
risks of large cascades and the incentives that shape them 
becomes ever more important. 

IV. MATERIALS AND METHODS 

A. Power grid topologies 

To understand coupling between multiple grids we 
requested data from the US Federal Energy Regulation 
Commission |49| . Using output of power simulations on 
many connected "areas" (distinct grids owned by differ- 
ent utilities) in the southeastern USA, we chose grid d by 
selecting the grid with the highest average internal degree 
and the grid, c, to which it had the most interconnections 
(8). Grids c, d have 439 and 504 buses (nodes) and aver- 
age internal degrees 2.40 and 2.91. For our purposes here, 
the only important details are the narrow degree distri- 
butions, the low clustering coefficients and the number of 
interconnections between the grids. Other details about 
the grids are in the SI. 

B. Toppling branch distributions 

To reduce the number of types in the branching pro- 
cess, we count the number of toppling events in each 
network rather than the number of grains of sand shed 
from one network to another. Here we derive the toppling 
branch distributions Ua , Ub from the shedding branch dis- 
tributions Qod- (For instance, Ua{ta, tb) is the chance that 
a toppled node in a causes ta,tb many nodes in a, 6 to 
topple at the next time step.) Note that a node topples 
if and only if it sheds at least one grain of sand. Thus 
a grain traveling from a network o to network d top- 
ples its destination node with probability 1 — (?o(i(0,0). 
Denoting (0, 0) by and the Binomial distribution by 
/3^(p) = (^)p'=(l -p)"-^ we have 

oo oo 

Ua{ta:tb) = Pai.ka,kb)PkS^ - 9aa (□)) A';" (1 - ^a, 

since the node must have at least ta many a-neighbors 
and at least tb many 6- neighbors, only ta,tb of which 
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topple (which are binomially distributed). The expres- 
sion for Ub is analogous. As an example, the probability 
generating function of Ua for Bernoulli-coupled random 
regular graphs R{za)-B{p)-R{zb) is 



Taking x = gives Sb, while taking x = 1 gives Sa- In 
practice, hundreds of terms of Sa,Sb can be computed 
exactly by truncating the sum in Q and using computer 
algebra systems. 



(P - PTa + {Za + l)(ra + " 1))"" (1 + piu " 1) + Zb) 



{Za + l)""- Za" {Zb + 1) 



C. Three methods for numerically solving the 
multidimensional self-consistency equations of a 
multitype branching process 

The most naive method to solve the self-consistency 
equations of a multitype branching process (such as Eqs. 
([3])) is to use a computer algebra system like Mathemat- 
ica or Maple to iterate ([s]) symbolically starting from 
Sa — Sb ~ I, expand the result, and collect coefficients. 
To obtain the coefficient Sa{ta,tb) exactly, it suffices to 
iterate ^ at least ta + h + l times. What is more, just 
a handful of iterations partially computes the tails (co- 
efficients of terms with high powers in Sa,Sb), but the 
amount of missing probability mass in the tails is unde- 
termined. This method was used in Fig. [3] 

A second method is to symbolically iterate Eqs. ([S]) 
at least ta + tb + 1 times and to use Cauchy's integral 
formula 
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where Z? C is a Cartesian product of circular contours 
centered at the origin, each of radius r smaller than the 
modulus of the pole of Sa closest to the origin [T3] . Then 
calculate one coefficient at a time using 
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A third method is Lagrange inversion, generalized to 
multiple dimensions by I. J. Good in 1960 |T3], a re- 
sult that has seen little use in the networks literature. 
We state the theorem in the language of the two-type 
branching process considered here, with the notation 
T = (rajTh), S = {Sa,Sb) and U = [UaMb]- For the 
slightly more general result that holds for arbitrary, fi- 
nite, initial population, see Ref. [13] . 

Theorem 1 [Good 1960] IfU{T) is analytic in a neigh- 
borhood of the origin andU{T) ^ (i.e., every type has a 
positive probability of being barren), then for x G {0, 1}, 
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D. Solving for the coefficients asymptotically (and 
why standard techniques fail) 

To determine the asymptotic behaviors of 
Sa{ta,tb),Sbita,tb) 'AS ta.tb OO, the trick is to 

solve for the inverses of Sa and Sb and to expand those 
inverses S~^ and Sj^^ at the singularities of Sa and iSf,. 
This technique works for isolated networks [THSl fT4] . 
but Sa and Sb have no finite singularities in for 
Bernoulli-coupled random regular graphs nor for power 
grids c, d. 

We demonstrate this failure of standard asymptotic 
techniques on the networks R{za)-B{p)-R{zb). Let a; :— 
S{f). Assuming an inverse S^^ of S exists, using ([3| we 
have 



where fi, v run over the types {a, b}, 5^ is the Kronecker 



delta, and 



is the determinant. 



■ = S-\u) = 



UaiL3) ' Ubiuj) 



(7) 



The generating functions Sa^Sb 
T*,fl*, respectively, if DS~^{t*) 

where the operator D — {d/duja,d/duJb 
partial derivatives. Differentiating Eq. [7\ 
the numerators to gives 



have singularities at 
= DS^'m = (0,0), 
is a vector of 
and equating 



on 



0, 



0, 



,dUa,^ 

-^^Ihb^^ 

on 



) = 0, 



(8a) 
(8b) 



The only solution to Eqs. ^ 

Zg+P-l P- Zb-l 
P - Zq - 1 ' p 

P- Zg-l zl+p-1 

p ' p- Zb-1 



is 



Ail = 



> (cx), 1 - Za), (9a) 
>{l-Zb,(x). (9b) 



This solution Q does not recover the singularities t* ~ 
(l,ci),/i* = (c2,l) (where ci,C2 are arbitrary) that we 
should obtain when the networks are isolated (p = 0) H . 
Moreover, although Eqs. (|8| vanish at the solutions (|9|, 
the derivatives of ([t]) do not vanish at these solutions 
(§, as DS-H^) = (oo,0), DS^\f2l) = (0,c3o). Thus 
we must discard solutions 

Solving only the left-hand equations in Eqs. Q yields 
singularities that do recover the correct singularities 
when p = 0: 



p-1 



{l-Za)ip- 



= Ma 



l-Za) 
p-1 



n 



{1 - Zb){p- I - Zb) J p-^o' 



^(l,n), (lOa) 



(aZ^,1), (lOb) 
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where Tf,, fj,a are arbitrary constants satisfying 

p- Za-l 



_ p - Zfc - 1 

Tb T : 

V 



P 



(11) 



SO that the derivatives of S^^{uj) ^ evaluated at ( 10 1 are 
finite. However, the derivatives on the right-hand sides 
of Eqs. (Isl) evaluated at the solutions ( 10 1 arc 



95-1 / p{zh + l){zl+p-l) [j§f^) 



(1 -Za)il-p + Za) {p {uJb - 1) + ^6 + 1)' 



and similarly for dS^^ /du}a{l^2) (interchange a and b). 
This derivative vanishes if and only if p = or p = 1 ~ -^a 
or Za = 1 or = oo. Hence we cannot find finite sin- 
gularities of Sa,Sb] these generating functions are en- 
tire functions with singularities only at Tb = fia — oo, 

Hence we have no singu- 



larities at which to do an asymptotic expansion for the 
coefRcients Sa,Sb, as one can for isolated networks 7- 

Other techniques exist for asymptotically approximat- 
ing the coefficients of generating functions, depending on 
the type of the singularity ([57 §5). Hayman's method 
f |57j §5.4) works for generating functions with no singu- 
larities in the finite plane (i.e., entire functions), like the 
Sa,Sb considered here. However, the theorem requires a 
closed form expression for the generating function, which 
we cannot obtain from the self-consistency equations for 
the synthetic and real interacting networks of interest. 
Developing techniques for asymptotically approximating 
the coefRcients of multidimensional generating functions 
with singularities at infinity, like those studied here, poses 
an important challenge for future studies of dynamical 
processes on multitype networks. 



V. APPENDIX: EFFECTIVE DEGREE 
DISTRIBUTIONS IN MULTITYPE NETWORKS 

Using the configuration model to generate multitype 
networks — including bipartite and multipartite graphs, 
graphs with arbitrary distributions of subgraphs, and the 
modular graphs considered here — requires matching edge 
stubs within and among types. For example, for the in- 
teracting networks considered here, the number of edge 
stubs pointing from network a to network b must equal 
the number from & to a. A standard practice in the lit- 
erature [IHl [SI] that is more restrictive than needed is 
to require that the averages of the inter-degrees over the 
degree distributions agree (e.g., (kab) — (kba) for two net- 
works a, b of equal size). In fact, most any degree distri- 
butions can be used, as long as conditioning on matching 
edge stubs among the types of nodes leaves some proba- 
bility. Requiring that the degree distributions satisfy, for 
example, (kab) = (kba): merely tips the scales in favor of 
valid degree sequences. 



Here we derive the effective degree distribution in mul- 
titype networks generated with the configuration model. 
The idea is simple: since nodes draw degrees indepen- 
dently, the probability distribution of the total number 
of edge stubs from, say, network a to 6 is given by a con- 
volution of the degree distribution. We state it for the 
two interacting networks considered here, but it can be 
easily generalized to, say, role distributions |54| . which 
require matching edge stubs with ratios different from 
one. 

Suppose networks a,b have Na,Nh many nodes, re- 
spectively. Let Kab, Kba be the random variables for 
the sequences of "inter-degrees" kab, kba of the nodes in 
networks a, 6, respectively, drawn from the input degree 
distributions Pa{kaa,kab),Pb{kba,kbb)- Suppose, for sim- 
plicity, that the internal and external degrees are inde- 
pendent {Pa{kaa,kab) = Paa{kaa)Pab{kab) , Pb{kba,kbb) = 

Pba{kba)pbb{kbb)) , although this is not essential. We de- 
note Y.k = h for fc e Z". 

Lemma 1 With the assumptions in the previous para- 
graph, the effective inter-degree distribution for network 
a is not the input one, Pr{Kab = k) = Yii'^iPabiki) , but 
rather the conditional probability distribution 



Pr{Kab = k I j:Kba = ^Kab) = Pr{Kab = k) 



Plai^k) 



where p*a^{-),pla{-) are Pab,Pba convolved Na, Nb times, 
respectively. 

Proof. Using the independence of Kab and Kba, we 
have 



PriKab = k I J:Kba = ^Kab) 



PrjKab = k) Pr(SXb„ = Sfc) 
Pr(Si?b, = ^Kab) 



In the numerator we recognize Y'i{Y,K}ja — Sfc) to be pba 
convolved Nb times, evaluated at E/c. In the denomina- 
tor, use independence, recognize convolutions and sum 
on ^ > 0. □ 
Lemma [l] shows that the effective inter-degree distri- 
bution of a is the input degree distribution pab reduced 
by an amount given by a fraction of convolutions. This 
reduction in probability governs how many invalid de- 
gree sequences must be generated before generating a 
valid one. For systems sizes on the order of 10* nodes, as 
considered here, generating degree sequences until pro- 
ducing a valid one is quite feasible, as it takes merely 
seconds. However, for millions of nodes or more, it is 
better to generate degree sequences Kab, Kba once, and 
then repeatedly choose a node uniformly at random to 
redraw its degree from its degree distribution until the 
degree sequences are valid. However, this method does 
not escape the effect in Lemma [T] which is often subtle 
but can be substantial if the supports of the convolutions 
of the degree distributions have little overlap. The inter- 
degree distributions used here, Bernoulli and Correlated- 
Bernoulli with identical expected total inter-degree, have 
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"much overlap" , so the effective inter-degree distribution 
is approximately the input one, and the correction factor 
in Lemma [T] can be neglected. 
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SI. 



Supporting information 



ROBUSTNESS OF OPTIMAL 
INTERCONNECTIVITY 



Sl.l. How p* depends on system size, connectivity, 
type of degree distribution 

Our central result is the minimum p* > in the chance 
of a large cascade in a network with a fraction p of its 
nodes connected to nodes in another network. Impor- 
tant questions are how p* depends on system size, be- 
haves in the thermodynamic limit (infinite system size), 
and depends on the internal degrees and on the degree 
distribution. Here we show that the qualitative form of 
Pr(Ta > C) in Fig. 4 appears to be generic. 

Doubling the system size and keeping the cutoff C and 
the dissipation / fixed (at C = half the number nodes in 
a and at / = 0.01) does not significantly change Fig. 4, 
because the dissipation / limits the chance of large cas- 
cades. But doubling the system size, doubling the cutoff 
C and halving the dissipation / (so that the largest cas- 
cades topple nearly the whole system) slightly decreases 
p* . Figure [ST] shows similar results as Fig. 4 but for a sys- 
tem half the size (1000 nodes/network, / = 0.02, cutoff 
C = 500). In the thermodynamic limit of infinite sys- 
tem size N , we expect p* to stay bounded away from 0, 
because what appears to be the determinant is the ra- 
tio of edges between the networks {pN) and within them 

{ZaN). 



pled random 4-regular graphs R{A)-B{p)-R{A)) increases 
their capacity and hence ability to withstand inflicted 
cascades, so the minimum in Y'i[Ta > C) is therefore 
larger, p* « 0.2 (Fig. S2), compared to p* « 0.12 for ran- 



dom 3-regular graphs with equally many nodes (Fig. SI I. 
The networks also have a wider range of optimal p, as ex- 
pected, given that pNa external edges is less significant 
than the ZaNa internal edges when Za = 4 (Fig. |S2| com- 
pared to Za = 3 (Figs. 4 and SI I. Nonetheless, the chance 
of a large cascade eventually increases with p. We also 
note that power grids more closely resemble random 3- 
regular than 4-regular graphs (see Table S2S2.1 below). 




Pr(r„„>500) 
Pr(r4„>500) 

Pi-(r„>500) 



FIG. S2: Increasing the internal degrees of both networks 
a,b from 3 to 4 (i.e., R(4)-B{p)-R(4)) slightly increases the 
minimum p* « 0.2 of the chance of a large cascade in network 
a. (Here: 10"^ nodes/network; / = 0.02; 10® grains dropped 
after 10^ transients.) 




1000 



Rank 



FIG. SI: For half the system size {Na = 1000), half the cutoff 
(C = 500) and double the dissipation (/ = 0.02) of the system 
in Fig. 4, the chance of a large cascade in a is qualitatively 
similar, and p* « 0.12 ± 0.02 is slightly larger (in simulations 
with 2 X lO'' grains dropped after 10^ transients). This gold 
curve and its critical point is stable to cutoffs C different from 
500 (200 < C < 800). Inset figures: rank-size plots on log- 
log scales of the largest cascades in network a (left) for p = 
10-^10-^10-^ and in power grid d connected to c by 0, 8 
or 16 edges. 



Increasing the internal degrees of nodes (say, to cou- 



Next we determined the effect of introducing some de- 
gree heterogeneity by adjusting interconnectivity p be- 
tween two Erdos-Renyi random graphs. Similar results, 
were obtained as for random regular graphs indicating 
that some degree heterogeneity does not affect these re- 
sults. (We did not test heavy-tailed degree distributions, 
since these rarely occur in the physical infrastructure net- 
works of interest, namely power grids.) 

Finally, we tested whether adding edges at random to a 
random regular graph produces an optimal connectivity 
that minimizes the chance of a large cascade. Specifi- 
cally, each node in a 3-random regular graph gains an 
extra, fourth edge stub independently with probability p 
(i.e., degree distribution P{k) is P(3) — l—p, P(4) = p). 
We find in simulations that Pr{T > 500) for a system 
of 10"^ nodes is largely independent of the mean degree 
z. In particular, the chance of a large cascade does not 
drop, as it does when initially adding edges between two 

This agrees with theoreti- 
the avalanche size 
We conclude that 



networks (Figs. 4, SI S2) 



cal results for isolated networks [7]: 
distribution s{t) ~ t^'^/'^ for large t. 
adding links between two networks affects vulnerability to 
large cascades in a different way than adding links within 
isolated networks. 
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SI. 2. Unstable p* for small cascades, stable p* for 
large cascades 

Networks seeking to mitigate cascades of small, inter- 
mediate or large sizes would seek different interconnec- 
tivity, as shown in the plots of Pr(£ < Ta < £ + 50) as a 
function of p in Fig. 6. In Fig. |S3| we show more plots to 
show the change in concavity at intermediate cascade size 
and the robustness of the location of the stable minimum 
p* « 0.075 for R{3)-B{p)-R{3). 



Pr(50 
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Pr(150 
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Pi(250 < T„ < 300) 
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FIG. S3: Plotsof Pr(^ < Ta <£ + 50), for £ = 50, 150, 850, 
as a function of interconnectivity p £ [0, 0.5], in simulations on 
R{3)-B{p)-R{3) with 2 x lO'' grains, 2 x 10^ nodes/network, 
/ = 0.01. Intermediate cascades (50 < Ta < 300) have an 
unstable critical point p* « 0.05, which changes concavity for 
cascades of size approximately 350. Large cascades Ta > 500 
have an stable critical point p* « 0.075. 



SI. 3. Increasing capacity fuels larger system-wide 
cascades 



Introducing new interconnections between networks 
slightly enlarges the largest global cascades. We ask to 
what extent is this due to the direction of links (pointing 
internally or externally) and due to the additional capac- 
ity for holding sand endowed by the new links (since ca- 
pacities of nodes are their degrees). To isolate these two 
effects, we run simulations on two random regular graphs 
interconnected by "Correlated-BernouUi" coupling: each 
node changes an internal edge stub to an external edge 
stub with probability p. 

The resulting total avalanche size distribution. 
Fig. |S1 S1.3[ does not significantly change for different 
values of interconnectivity p. Thus, what causes the 
slight enlargement of the largest global avalanches is the 
slight increase in total capacity of the network, not the di- 
rection of links. Moreover, the enlargement of the largest 
global cascades is consistent with the increase in capac- 
ity. For the simulations on random 3-regular graphs in 
Fig. |S1 S1.3[ the nth largest avalanche contains on aver- 
age 1.4±2.6 (mean ± standard deviation) more topplings 



with p = .1 compared to p = .005, which is an insignifi- 
cant difference. By contrast, in an analogous simulation 
on i?(3)-i?(0.1)-i?(3), in which we introduce extra hnks 
with probability p, the nth largest avalanche contains 
204.8 ±92.5 more topplings than i?(3)-S(0.005)-i?(3), a 
significant difference for networks with 10"^ nodes. Fur- 
thermore, this difference is on the order of the additional 
capacity of 2 x 10^ x (0.1 - 0.005) = 190 among 2 x 10^ 
nodes with B{.1) compared to B{.005) coupling. 
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FIG. S4: When nodes change exactly one of their inter- 
nal edge stubs to be an external edge stub with probabil- 
ity p (Regular(2:a)-Correlated-Bernoulli(p)-Regular(2;6) rather 
than receiving an external edge stub with probability p 
(Regular(2a)-Bernoulli(p)-Regular(2i,)), global cascades are 
not significantly amplified. This indicates that it is the small 
increase in capacity (due to the additional edges and the fact 
that capacities of nodes are their degrees) rather than the di- 
rection of edges (internal versus external) that causes the in- 
crease in global cascades with increased Bernoulli-distributed 
coupling between networks. The main plot is the total 
avalanche size distribution s{t) in a simulation with 10® grains 
of sand dropped (after 10^ grains dropped without collect- 
ing statistics, begun from initial amounts of sand chosen uni- 
formly at random from zero to one less a node's degree) on 
two networks with 10^ nodes each, with dissipation of sand 
/ — 0.05 and Za = Zb = 3. Inset: rank-size plot of largest 
5 X 10'' cascades, which are nearly indistinguishable. 



We find similar results when introducing eight addi- 
tional interconnections between power grids c and d. The 
nth largest global cascades in simulations are larger by 
an amount on the order of the additional capacity of the 
two networks. Thus the additional capacity due to the 
new interconnections, not the direction of links, largely 
explains the amplification of system-wide cascades. 



SI. 4. Bounding the critical load disparity r* 
power grids c, d 



for 



In the main paper, we showed that for 16 intercon- 
nections between power grids c and d, the critical load 
disparity is r* sa 15. That is, when sand is dropped 15 
times more frequently on c-nodes than on d-nodes, the 
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largest inflicted cascades from c to d approximately equal 
in size the largest local cascades begun in d. (Recall that 
by "inflicted cascade from c to d" we mean the number of 
topplings in d in an avalanche begun in c, and by "local 
cascade in d" we mean the number of topplings in d in an 
avalanche begun in d.) If d builds more interconnections 
(or increases the load disparity r), the largest inflicted 
cascades become larger than the largest local cascades. 
On the other hand, delete interconnections (or decrease 
the load disparity r), and d could mitigate its largest 
local cascades more than the enlargement of the largest 
inflicted cascades. Hence r* — 15 and 16 interconnec- 
tions represent a "modularity equilibrium" in local and 
inflicted cascades. 

Here we use simulations to approximately bound the 
critical load disparity 10 ^ r* ^ 20 for 16 interconnec- 
tions between c and d. As shown in Fig. SI S1.4[ for 
r = 10, 20 the largest inflicted cascades (green dashed) 
are smaller and larger, respectively, than the largest lo- 
cal cascades (gold dashed), indicating that 10 ^ r* ^ 20 
for 16 interconnections. For 8 interconnections (solid 
blue and red curves), the critical load disparity r* is ev- 
idently larger than 20, as the largest inflicted cascades 
are smaller than the largest local cascades for r = 20 
(red curve below blue curve). 
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FIG. S5: Same plot as in Fig. 7, but for load disparity r = 10 
(main plot) and r = 20 (inset). For 16 edges between grids 
c and d (dashed curves), the largest inflicted cascades from 
grid c to d (green dashed curve) are slightly smaller than the 
largest local cascades from d to itself (gold dashed curve) for 
r = 10, indicating that r* ^ 10. For r — 20 (inset), the 
largest inflicted cascades are slightly larger than the largest 
local cascades, indicating that r* ^ 20. 



not appear to contain multiple electrical grids connected 
to one another. As a result, we requested data on power 
grid connectivity from the US Federal Energy Regulation 
Commission (FERC) via the Critical Energy Infrastruc- 
ture Information program [3]. We focused on the south- 
eastern region of the United States, for which wc had 
the output files of power simulations of various grids. 
The southeastern region consists of areas, distinct grids 
owned by different utility companies. The largest areas 
contain thousands of buses, an electrical grid term for the 
connection points that link generators, transmission lines 
and transformers. We ignore wind turbines, for they do 
not belong to specific areas. Among the areas with at 
least 100 nodes, we chose grids c and d by first selecting 
the grid with the highest average internal degree (2.91), 
and then selecting the area, c, to which it had the most 
interconnections (8 of them) . Grid d consists of one giant 
component with 504 nodes, while grid c consists of one 
giant component of size 439, and we ignore 14 additional 
nodes. 

Two of the statistics of power grids c and d are im- 
portant to our study: the narrow degree distribution 
and small clustering coefficient. Some details are in Ta- 
ble S2S2.1 The average clustering coefficient (C), the 



fraction of possible triangles that exist, begins to mea- 
sure how locally tree-like power grids are. We find that 
(C) « 0.05 is low, yet an order of magnitude larger than 
the clustering coefficient of a random 3-regular graph 
with as many nodes as grid d. The average shortest path 
length among pairs of nodes, {£), which has been iden- 
tified as a likely source of difficulty for predicting cas- 
cades on networks |5], is rather large in the power grids 
((£) « 10), due to their two-dimensional, spatial, nearly 
planar structure. By contrast, random 3-regular graphs 
have significantly smaller diameter {£) 7.09 ± 0.03. 





c 


d 


c & d 


3- Regular 


# nodes 


439 


504 


943 


504 


# internal edges 


527 


734 


1261 


756 


# external edges 


8 


8 






(^internal) 


2.40 


2.91 


2.69 


3 


(^external) 


.0205 


.0179 






(C) 


0.0109 


0.0821 


0.0488 


0.003(2) 


{() 


9.32 


8.26 


11.42 


7.09(3) 



TABLE SI: Summary statistics of power grids c and d (in 
isolation and coupled together) and of a random 3-regular 
graph ( "3- Regular" ) with the same number of nodes as d. The 
statistics for the random 3-regular graph are averages over 
1,000 realizations, with the standard deviation in parentheses 
to convey the fluctuation in the last digit. 



S2. POWER GRID TOPOLOGIES 



S2.1. Power grids c and d 



S2.2. Correlation between degree and load 



The readily available datasets on power grids — IEEE 
test cases P, Western States |2, and NYISO 0— did 



In the sandpile model studied here, we choose the ca- 
pacity of every node to be its degree, so that toppled 
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nodes unambiguously shed one grain to each neighbor. 
To examine the reasonableness of this assumption for real 
infrastructure, we sum the apparent power on the edges 
incident to each bus in power grids c and d, and we plot in 



Fig. S2 S2.2 a density histogram of apparent power versus 
degree of the nodes. Most buses have low degree and low 
apparent power, and some buses have high degree but low 
power, while others bear high power among few connec- 
tions. But the general trend is positive: each additional 
degree correlates with 123 MVA of additional apparent 
power, though the correlation is tenuous (i?^ = 0.30). 
This suggests that using the sandpile model with capac- 
ity equal to degree in order to study cascades of load in 
power grids — and by extension other infrastructure — is 
not unreasonable for the purpose of obtaining heuristic 
understanding. 
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FIG. S6: Apparent power versus node degree for power grids 
c and d, showing a weak correlation between the degree of 
nodes and their load. The best fit line power = —17 A -f 
123. 8fc indicates that each additional degree k correlates with 
an increase in apparent power of 123.8 MVA [B^ = 0.30). 



S3. BRANCHING PROCESS APPROXIMATION 



S3.1. Comparing theory and experiment 



we show results for the same branching process predic- 
tion (seven iterations of the self-consistency equations, 
started from Sa ^ = 1) with two independent sim- 
ulations with dissipation rate of sand / = 0.05. Halv- 
ing the dissipation rate noticeably extends the tails of 
the marginalized avalanche size distributions — i.e., the 
largest avalanches become larger. Iterating the self- 
consistency equations an eighth time would more accu- 
rately compute the probability mass in the tails, but is 
beyond our computer resources. Showing the results of 
two independent simulations (blue and green curves of 
Fig. S3 S3. 1 1 illustrates the variation in the cascade size 
distributions. In particular, variation between simula- 
tions is apparent only in the tails of the distributions, as 
the largest avalanches are so rare. Even though the be- 
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FIG. S7: Blue and green curves: independent simulations of 
10'^ grains of sand dropped on power grids c and d, labeled 
a and b respectively to match our theoretical notation. As 
in Fig. 3, we collect statistics after 10^ grains are dropped, 
initialized with amounts of sand chosen uniformly at random 
from zero to one less a node's degree, with dissipation of sand 
/ — 0.05 (half the value used in Fig. 3). Red curve: branch- 
ing process approximation using the empirical degree distri- 
butions of c and d, obtained by symbolically iterating the 
self-consistency equations (Eq. (|3|) seven times, expanding 
the result and collecting coefficients. The tails are decreased 
in the simulations due to finite system size (439 in c, 504 in 
d) and the dissipation of sand /, while the tails of the branch- 
ing process approximation miss probability mass due of only 
iterating seven times (which took a week on a 2GHz laptop) . 



We choose the dissipation of sand / (the chance that 
a grain of sand is deleted as it is shed from one node to 
another) so that the largest cascades in simulations top- 
ple almost all the network. (As a rule of thumb, take 
/ = 20/N, where N is the total number of nodes.) De- 
creasing or increasing / extends or shortens the tail of the 
avalanche size distribution, respectively. Thus we tune / 
to utilize the entire system and to achieve a power law 
over as many scales as possible, since we see power law 
behavior in, e.g., blackouts in power grids [B]. 

In Fig. 3 of the main paper we compare simulations of 
sandpile cascades with theoretical predictions for dissipa- 
tion rate of sand / — 0.1. This somewhat large dissipa- 



tion rate mitigates the largest cascades. In Fig. S3S3.1 



havior of sandpile cascades on single networks is rather 
robust to network structure [7Hl2). we see sensitivity to 
the sparse connections among modules, as demonstrated 
in the main paper. What is more, the particular de- 
grees of nodes with external links profoundly affect the 
inflicted avalanche size distributions, J2t ^aita^h) and 
X^tj *f)(^a' ^b)- Notice the characteristic "blips" (devi- 
ations from straight lines) that occur in the blue and 
green curves (the simulations) in the inflicted avalanche 
size distributions in the top-right and bottom-left plots 
in Fig. |S3 S3.1| For instance, there is a particularly high 
probability that an avalanche begun in d topples 8, 9 or 
10 nodes in c. We suspect this is due to the outlier node 
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in network c that has one external link and ten inter- 
nal links — far more internal links than all other c-nodes 
with an external link. Consequently, a cascade in d that 
leaks across to c via this node topples with high proba- 
bility around 8, 9 or 10 nodes in c because of this node's 
high internal degree. Similarly, the blip at size 3 of the 
distribution of inflicted cascade size from c to c? results 
from the average internal degree 2.5 of d- nodes with an 
external link. It is surprising, given the indifference of 
sandpile cascades to network structure, that the particu- 
lar internal degrees of nodes with an external link greatly 
affect the marginalized avalanche size distributions. 
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FIG. S8: Comparison of theory (red) and simulation (blue) of 
the chance of toppling 10 nodes in a and < a; < 30 nodes in 
b, in log-log scales. The simulations were Bernoulli-coupled 
random 3-regular graphs with lO"' nodes each, dissipation / = 
0.02, and 2 x 10® grains dropped. 



S3. 2. Distance between the avalanche size 
distributions and the product of their marginals 



S3. 4. Capacity disparity 



When we compare the branching process predictions 
with simulations, such as in Figs. 3 and |S3 S3.1[ we re- 
duce the dimensions of the joint Sa(ta,tii) by marginal- 
izing in order to plot one-dimensional curves. This 
raises the question: how far are the products of the 
marginals, such as Sa(ta, tb))(I]t^ Sa(ia, 4)), 

the joint distributions? Using the first 11^ coefficients 
of Sa{ta,tb) (for < ta,tb < 11) computcd using La- 
grange inversion (Theorem [l]), we computed the distance 
between Sa{ta,tb) and the product of its marginals — in 
the 2-norm, 1-norm and KuUback-Leibler divergence — as 
a function of Bernoulli coupling p between two random 
3-regular graphs. With increasing coupling p, the cas- 
cade sizes ta,tb become increasingly correlated, and the 
joint distribution Sa grows increasingly distant from the 
product of its marginals. This suggests that comparing 
marginals, as in Figs. 3 and |S3 S3.1[ is an insufficient test 
of the branching process, and thus we next subject the 
branching process to its most stringent test. 



S3.3. 



Fine-grained comparison of branching 
process and simulation 



Since marginalizing Sa,Sb, as in Figs. 3 and |S3 S3.1[ 
may obscure deviation between theory and the branching 
process, we next compare the joint distribution Sa{ta, tb) 
to simulation. In Fig. |S3 S3.3| we plot Pr(Ta = 10, Tb = x) 
for < a; < 30 in the simulation (blue) and branch- 
ing process prediction (red), for interconnectivity p = 
0.005,0.01,0.1. We compute the branching process pre- 
diction by calculating Sa{ta,tb), Sb{ta,tb) for ta = 10 and 
< tfe < 30 using Lagrange inversion (Theorem [I]). Since 
sand is dropped on nodes uniformly at random and each 
network has 10^ nodes, the avalanches begin in a or in 
b with equal probability, so that Pr(Ta — 10, Tf, = x) = 
(sQ(10,a;) -|- Sb{l0,x))/2. Although the branching pro- 
cess under-predicts theory for larger avalanches (espe- 
cially with greater interconnectivity) , we see qualitatively 
good agreement. 



Cascades inflicted from large capacity networks to 
small capacity networks are larger than those from small 
to large capacity networks. Here we derive the heuris- 
tic formula for this effect (Eq. ([4|) using the multitype 
branching process, and we show qualitative agreement 
with simulation. 

Differentiating the self-consistency equations ([s]) with 
respect to and Tb and setting = t;, = 1 yields four 
equations for the first moments of the avalanche size dis- 
tributions Sa,Sb' 
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(S12a) 
(S12b) 



with common denominator 



7 = -1 + Mb Ma - Ma (Mb - 1) + Mb ■ (S13) 

For the first moments of the branch distributions Ua , u& 
for R{za)-Bip)-R{zb), 
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l + Zb 



(S14) 



(and vice versa for Ub, interchanging a with 6), the de- 
nominator (S13) is zero, and Eqs. (S12) are all infinite. 



This is not a surprise, because in isolation the networks' 
cascade size distributions are asymptotically power laws 
with exponent —3/2, which have infinite first moments. 
However, one can compare the rates of divergence of the 
mean infiicted cascade sizes, Eqs. (S12b|, by computing 
then ratio, {sa}b/{sb)a = (1 + Za)/(1 + Zb). No other 
ratios appear to contain useful information (nor depen- 
dence on p with which to estimate the critical amount of 
interconnectivity p* ) : 



\^a)a \^a)a 

{Sb)a {Sb)b 



Za 
Zb '' 



{Sa)a 



= 1. 



In all simulations on R{za)-B{p)-R{zb) with Za ^ Zb, 
the heuristic formula Eq. (|4| works qualitatively, in that 
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inflicted cascades are always larger from the larger ca- 
pacity network (i.e., Za < Zb =^ {sa)b < {sb)a)- 
But finding quantitative agreement with Eq. Q is dif- 
ficult because the avalanche size distributions are power 
laws with such heavy tails that their variances diverge. 
Figure [SSSO] is an example: for i?(3)-B(0.005)-i?(3), 
{sa)b = 0.086 ± 0.019, and {sb)a = 1-32 ± 0.023 (mean ± 
standard error), which yields {sa)b/ {sb)a = 0.65 ± 0.019, 
while theory predicts {1 + Za) / {1 + Zb) =4/5 (20% error). 
Nevertheless, Eq. Q provides a simple, useful heuristic 
for the dangers of small capacity networks connecting to 
large capacity networks. 



tiating expressions raised to large powers, so the largest 
coefficients are the most difficult to compute. 

The Cauchy formula, Eq. ^ of the Materials and 
Methods, uses integration rather than differentiation, 
which can be more numerically stable [14j . But com- 
puting the largest coefficients Sa{ta,tb) with this method 
requires integrating an increasingly large expression, 
namely the result of iterating the self-consistency equa- 
tions at least ta+tb + l times, starting from Sa ^ Sb = I- 
This, too, takes a long amount of computation time: 
computing Sa{ta,tb) for < ta,tb < 20 on a typical lap- 
top computer would take on the order of a week. 
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FIG. S9: Qualitative agreement between the theoretical 
prediction (Eq. Q) and simulation on i?(3)-B(0.005)--R(4) 
(2 X 10® grains, 10^ nodes/network, / = 0.02). Main plot: 
marginalized avalanche size distributions of inflicted cascades 
from o to 6 (blue), 6 to a (red). Inset: rank-size plot of the 
largest 10'' inflicted cascades from a to fe (blue), 6 to a (red). 




S3. 5. Computational challenges in showing the 
locally stabilizing effect 

Although the locally stabilizing effect of interconnec- 
tions is apparent in simulations, demonstrating it an- 
alytically poses mathematical and computational chal- 
lenges. Here we use multidimensional Lagrange inver- 
sion (Theorem [l] in the Materials and Methods) to solve 
for the probabilities of the smallest avalanches, Sa{ta,tb) 
for < ta,tb < 10. Next we plot in Fig. IS3 S3.5| A the 
marginalized avalanche size distribution X]tb=o Sa{ta,tb). 

Although it appears in Fig. |S3 S3.5| A that the largest 
avalanches become less likely with increasing intercon- 
nectivity p, in accordance with the results in simulations, 
this figure has a caveat: we have only computed up to 
avalanche size at most 10 in each network. As a result, 
the right-hand tail of Fig. |S3 S3.5| ^ lacks probability mass 
compared to the actual marginalized avalanche size dis- 
tribution J2ti!=Q Sa(ta,tb), because, for example, there is 
a significant chance of an avalanche of size = 10 and 
tb = 11, 12 or 13. This figure could be improved by com- 
puting Sa{ta,tb) along "vertical strips" < tb < t^°-^ , 
where ^ 1, for just a handful of values ta- However, 
computing the large coefficients of Sa requires differen- 



FIG. SIO: Panel A: Increasing the interconnectivity p between 
two random 3-regular graphs appears to mitigate the tails of 
the local avalanche size distribution "^-/.^ Sa{ta,tb)- Plotted 
are the first 11 x 11 coefficients of Sa computed using multi- 
dimensional Lagrange inversion [TJ] for p = 0,0.5, 1. Panel 
B: Increasing the interconnectivity p between two random 
3-regular graphs "smears out" the avalanche size distribu- 
tion Sa{ta,tb), SO that cascades large in both networks be- 
come more likely. Plotted are Sa{ta,tb) for < ta,tb < 10, 
computed using multidimensional Lagrange inversion [13' for 
p = 0,0.1, 1. We plot zero probabilities in white and positive 
probabilities in a logarithmic color scale from orange (high 
probability) to bright red (low probability). 



Nonetheless, even the smallest 11^ coefficients of Sa 
and Sb are useful. Fig. |S3 S3.5p , for example, shows the 
joint avalanche size distribution Sa{ta, tb) for two random 
regular graphs with Bernoulli(p)-distributed coupling, for 
p — 0,0.1,1. As the interconnectivity p increases, the 
cascades begun in a increasingly "smear out" among a 
and 6, and cascades are more frequently large in both 
networks at once. 
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